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Abstract 



We present a new simulation scheme based on the Lattice-Boltzmann method 
to simulate the dynamics of charged colloids in an electrolyte. In our model 
we describe the electrostatics on the level of a Poisson-Boltzmann equation 
and the hydrodynamics of the fluid by the linearized Navier-Stokes equations. 
We verify our simulation scheme by means of a Chapman-Enskog expansion. 
Our method is applied to the calculation of the reduced sedimentation velocity 
U /Uq for a cubic array of charged spheres in an electrolyte. We show that 
we recover the analytical solution first derived by Booth (F. Booth, J. Chem. 
Phys. 22, 1956 (1954)) for a weakly charged, isolated sphere in an unbounded 
electrolyte. The present method makes it possible to go beyond the Booth 
theory, and we discuss the dependence of the sedimentation velocity on the 
charge of the spheres. Finally we compare our results to experimental data. 
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I. INTRODUCTION 



The simulation of the dynamics of colloidal suspensions is a challenging task. The rea- 
son is that the movement of the colloidal particles can be on a time scale which is orders 
of magnitudes slower than that of the solvent particles (e.g. seconds versus picoseconds). 
Therefore, simulation methods such as Molecular Dynamics, that account for the fluid par- 
ticles explicitly, are not well suited to study the dynamics of colloidal suspensions because 
one would spend most of the simulation time solving the equations of motion of the fluid 
particles. The situation for charged colloids is even worse because one has to take into 
account the long-ranged Coulomb interactions which are already very time consuming in 
simple ionic liquids. 

One possibility to circumvent these problems is to avoid explicit simulation of the sol- 
vent particles and describing the interaction between the colloidal particles by means of an 
effective potential. In the case of charged colloidal systems this is usually a Yukawa-like 
potential which gives in many cases quite an accurate description of interaction-dependent 
properties |jl|J^. But the approach with effective interactions neglects completely the hy- 
drodynamic interactions between the colloidal particles which stem from the momentum 
transport through the solvent. In order to take this into account one has to treat the hy- 
drodynamics of the solvent at least on a coarse-grained level. An effective scheme, that was 
developed to solve efficiently the Navier-Stokes equations, is the so-called lattice Boltzmann 
method (LBM). The LBM is a pre-averaged version of a lattice g Boltzmann equa- 

tion is solved on a lattice such that the Navier-Stokes equations are recovered (reviews of 
the method are the Refs. Recently, the LBM was applied to simulate the dynamics 

of colloidal systems, such as the rotational and translational short-time dynamics of colloids 
||6HlO|, the diffusion of colloidal particles in confined geometry and the dynamics 

in porous media [|T^-|T^. Also other complex systems like polymer solutions [0 have been 
investigated by LBM. 

Several attempts to apply the LBM to charged systems have been reported in the lit- 
erature. He and Li [19| proposed a scheme which is appropriate to study electrochemical 
processes in an electrolyte. However, in this method it is assumed that the fluid is locally 
electrically neutral which cannot be true for the part of an electrolyte forming the electrical 
double layer around a macroion. Thus, the LBM of He and Li cannot be used to describe 
the dynamics of suspensions of charged colloids. A different LBM for charged systems was 
suggested by Warren [^]. The central idea of this method is the introduction of external 
charge densities ps for the ionic species of type s which are propagated with the one particle 
distribution function of the LB equation by means of the so-called moment propagation 
method These ionic densities are coupled back to the mass current of the LB equation 
via a chemical potential which consists of a term oc Inp^ and a term proportional to the 
electrostatic potential such that it only gives a contribution if the ion densities are not dis- 
tributed as expected in equilibrium by a Boltzmann distribution (e.g. because of an external 
electrical field). Finally, the electrostatic potential is determined from the charge densities 
by means of a Poisson equation solver. The main drawback of Warren's method is that the 
charge densities are introduced as additional physical quantities which are independent from 
the mass density in the LB equation, and thus, this method is not self-consistent and the 
correctness of coupling of the charge density to the mass current is not obvious. 
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Our LBM for charged colloidal suspensions is inspired by Warren's approach, but it does 
not introduce ionic species as additional quantities. The details of our method can be found 
in the following sections. We apply it to the determination of the sedimentation velocity of 
an array of charged spheres in an electrolyte. 

The article is organized as follows: In Sec. II we give a brief introduction into the 
electrokinetic equations of motion with which we model the dynamics of the fluid. In Sec. Ill 
the LBM is presented which solves the latter equations of motion on a lattice. The LBM is 
verified in Sec. IV by means of a Chapman-Enskog expansion. And in Sec. V we show our 
results for the sedimentation phenomena. We conclude with a discussion of the method and 
the results. 



II. THE ELECTROKINETIC EQUATIONS 

In this section we introduce the equations of motion for a hydrodynamic description of 
charged colloidal suspensions. These equations can be found in standard textbooks (see 
e.g. Ref. |22|). 

Consider a system of macroions with radius a in an electrolyte consisting of two ionic 
species which have charges +Zie and —Z2e, respectively (where zi and Z2 are the valencies 
of the ions and e is the charge of a proton). The densities of the ions, pi{r,t) and p2{r,t), 
are conserved quantities and therefore each of them follows a continuity equation: 

^ = -V-J. . = 1,2. (1) 



The current Jg is given by 



Js = PsU- DsVps - ZsDsPsV^ . (2) 



The first term in (^) in which u denotes the flow velocity is the convection current whereas 
the two other terms describe the diffusive current and the current due to the electrostatic 
potential $. Dg denotes the diffusivity of ions of type s and $ is the electrostatic potential 
in dimensionless form, $ = -j^^ {ks'- Boltzmann constant, T: temperature). 

$ is determined by the Poisson equation. 



V^^ = -ATxlB\^^z,p, + a^ (3) 
where the Bjerrum length Ib is defined by 

'--sir- W 

e is the dielectric constant. At a distance Ib the Coulomb energy of one ion due to another 
ion is equal to ksT. a denotes the charge density of the macroion. We assume that the 
charge Z of each macroion sits on its surface in the form of uniformly distributed point 
charges. 
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If we set M = in Eqs. (|l]) and (0) we yield the equilibrium solution for the ion densities 

as 



Ps{r) = Ps exp (-z,<l>(r)) • (5) 

By putting the Boltzmann distribution (^) into the Poisson equation this leads to the 
so-called Poisson-Boltzmann equation in which correlations between the ions are neglected. 
Moreover, if one linearizes the exponential function in Eq. (|^) (Debye-Hiickel theory) it is 
possible to solve Eq. @ analytically and the result is a Yukawa potential, 

= K'-^^t^^ . (6) 



|r| 



The so-called Debye length is defined by 



^ = ^ . (7) 



The potential (j^) is even a good approximation in the non-linear case but the prefactor K 
is then different from the one of the linear theory. 

The equation of motion which has to be specified finally is that for the total mass current 
of the fluid, pu = (Es=i Ps + Pn)u, where p„ is the density of the neutral part of the fluid. 
We assume that our fluid can be described by the linearized Navier-Stokes equations for low 
Reynolds number flow. Hence, the equation for pu with a body force due to the electrostatic 
potential is 

^ = uV'pu - Vp - ksT J2 ^.Ps V$ (8) 

where p is the pressure and u is the kinematic viscosity. If the equation for the pressure is the 
one of an ideal gas, p = ksTp, p can be decomposed into an electrostatic and a neutral part, 
Pe = ksTY^I^^ ps and Pn = ksTpn, respectively. The sum of — Vpe and the electrostatic 
body force, i.e. the last term in Eq. (|^), is zero, if the ion densities have relaxed to their 
equilibrium distribution, Eq. (H). 



III. THE LATTICE BOLTZMANN METHOD FOR CHARGED COLLOIDS 

We have developed a simulation method to solve the non-linear coupled Eqs. (P, @, 
and (|^) numerically. For this we use concepts which are well-known from the LBM. 

In this method the discretized version of a Boltzmann equation is solved numerically on 
a lattice on which every lattice point represents a cell of particles. The central quantity is 
the one-particle distribution function ni{f,t) which describes the number of particles on a 
lattice node r at time t with a discrete velocity Cj. The discrete space of velocities {cj} is 
chosen such that no artificial anisotropic terms appear in the corresponding equations in 
the continuous limit. In our case the velocity space consists of 18 vectors of which from 
a given lattice node 6 point to the nearest and 12 to next-nearest neighbors on a simple 
cubic lattice. This velocity space can be constructed by projecting the unit vectors of a 
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four dimensional FCHC lattice onto three dimensions. It is one possible choice of a velocity 
space which exhibits the required isotropy. 

The equation of motion for nj(r, t) consists of two steps, a collision and a propagation 
step. In the collision step the interaction between the particles is taken into account which 
results in the post-collision function n*{f,t*) at the collision time t*. In the propagation 
step riiif.t) is updated by 

ni{r + Cut + I) = n*{r,t*) . (9) 

In this equation the lattice constant, the time step, and the mass of a particle is set to unity. 
The density p(r*, t) and the mass current j = pu are given by the zeroth and first moment 
of ni, respectively, p = Ei ni{f, t), j = J2i ni{f, t)ci. 

In the case of the charged system a one-particle distribution function for each ion species 
and a neutral part is required. The purpose of the neutral part is to keep the viscosity 
essentially constant through the fluid. Thus, it is chosen such that its value at a given 
lattice point is much higher than that of the ionic densities. We make the following Ansatz 
for the post-collision function n^* for the counter- and colons, s = +, — , respectively, and 
the neutral part, s = n (we also take into account rest particles by the index i = 0): 

-r(r, n = ^Ps{r, t) [l + ^J^^M t) ■ q) , (10) 

<(f,t^) = (l-7,)p,(f,t) . (11) 

The factor Wj is a weighting factor which is equal to 2 for the q in the direction of nearest 
neighbors and equal to 1 for the remaining q. So it satisfies the normalization constraint 
J2iti U = 1- the following we define also wq := 0. With the parameter < 7s < 1 the 
diffusivity Dg of the particles of type s can be varied. The latter quantity is given by 

Ds = f 7s (12) 

which is shown in the next section. Csv is the sound velocity which is 1/ \/2 for our model . 
The density p' is defined by p' = J2s1sPs- 

The propagation step for our charged system is 

n|(f+Q,t + l)=<(r,0, (13) 
n^(f,t + l)=<(f,0 . (14) 

Different propagation rules have to be established at the surface of the macroions and at 
walls. Here we use the bounce back rules suggested by Ladd 0] which lead to no-slip 
boundary conditions. In this scheme one puts a sphere which represents a macroion onto 
the lattice whereby its surface cuts links between lattice nodes. The boundary nodes are 
defined halfway along these links and the population functions nf which point to the direction 
of the boundary nodes are reflected back during the propagation step. In the case of moving 
boundaries there is a momentum transfer between the boundary nodes and the fluid. In this 
case the self consistent scheme derived by Lowe et al. ||^,|T^ can be used. 

The aforementioned way of mapping a sphere onto the lattice introduces fluid inside and 
outside the sphere. In our scheme we assign the charge of the macroions in that the inner 
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fluid is an electrolyte with net charge Z. Charge neutrality requires that the total charge 
of the outer fluid equals the one of the inner fluid of the macroions. Of course, it is not 
allowed in our scheme that outer fluid leaks through the surface of the sphere. Therefore, 
only small movements of a macroion are possible such that the center of mass of the sphere 
can be fixed, and only a momentum transfer with the fluid takes place. 

The densities ps and the total mass current j cannot be inferred simply from the zeroth 
and first moments of the ra^'s because we have to take into account their coupling to the 
gradient of the electrostatic potential. If ps and j are calculated as follows, 

18 

p,(f,t + 1) = ^ (n|(f,t + 1) - |^p,(f - cl,t)Vl'(f - c-;,t) ■ cl) (15) 



and 




j(r, t + 1) = E E t + - cizslsPsir, t + l)V$(f, t + 1) . (16) 



we are consistent with Eqs. (|l|) and (||) in the continuous limit. does not couple to the 
neutral part of the fluid. This is guaranteed in Eqs. ( PTBD and (0) by setting Zs = for 
s = n. 

If one replaces nf{r, t+1) by n^{r—Ci, t) in Eqs. ([l5|) and ([T6|) by using Eq. ( [l3|) it becomes 



clear that Eqs. (|T5D and (|T^) are the discrete versions of Eqs. (|l]) and (H), respectively. 
Thereby, the second terms in Eqs. (|l^) and ([iBD correspond respectively to the current due 
to the electrostatic potential in Eq. (|1|) and the electrostatic body force in Eq. (H). We show 
this explicitly in the next section by means of a Chapman-Enskog expansion. Of course, 
with the V$ terms in Eqs. (0) and (|TB|) the conservation of the ionic densities and the mass 
current is still fulfilled. This is guaranteed because $ is determined self-consistently from 
the ionic densities by means of the Poisson equation. 

This means that we have to solve the Poisson equation at each time step in order to 
determine the electrostatic potential from the ionic densities p+ and p_. For this purpose 
we use a successive over-relaxation (SOR) scheme in which one looks in principle for the 
stationary solution of a diffusion equation |j23|. But in contrast to a normal diffusion equation 
one introduces an acceleration parameter 1 < uo < 2 such that the potential is obtained as 
the iterative solution of the following equation: 



+ {1 - u;)^h{r,t) . (17) 



In this equation we have denoted the iteration time by h (the unit for an iteration step is 
again set to unity). For a sufficient number of iteration steps converges to $. We have 
found that a suitable choice for u is 1.45 which guarantees stability and optimal acceleration. 
The gradient of $ which we need in Eqs. (|15]) and ([T6|) is given by 

v<i(f) = -y: ^h^- m ■ (18) 



Note that we use in Eqs. ([T^) and (|T8[) the same discrete space as the velocity space in the 



Lattice-Boltzmann equations. This means that the truncation error which is caused due 
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to the discrete representation of the derivatives is of fourth order in contrast to a second 
order truncation error on a simple cubic lattice with 6 vectors {cj} pointing to the nearest 
neighbors from a given lattice node. The stability and efficiency of the SOR algorithm are 
further optimized by making use of a partially decoupled red-black Gauss-Seidel scheme 

El. 



IV. CHAPMAN ENSKOG EXPANSION 

By means of a Chapman-Enskog expansion we show in this section that the set of 
discrete equations from the previous section indeed recover the Eqs. (||), (0), and (P) in the 
continuous limit. 

As a first step we rewrite Eqs. (|T5|) and ( pJf ) for the densities Ps and the partial currents 

Is = IsPsJIp'- 

Ps{r,t) = Zl"^^ (ps(f-Q,t- 1) +Js{f-Ci,t- 1) -Ci 

+^Ps{r-Ci,t-l)VHf-Ci,t-l)-c,^ , (19) 
Js{r,t) = {ps{r- Q,t- l)ci + js{r- Cut- 1) ■ QCi) 

i 

-cl^,ZsPs{f,t)V^f,t) . (20) 

If we now expand the functions of the form f{r — Ci, t — 1) (/ = Ps,js, PsVl>) around position 
f and time t up to second order, 

fir- Ci, t - 1) = /(r-, t) + (^-dt - Q,V, + + Q„Va)') /(r, t) , (21) 

we obtain the following equations for ps and js, 
1 

dtPs = -dlps -V -l + dtV -1 + (vVs + ZsV ■ PsV^ + ZsdtV ■ p,V<l) (22) 

dtl = -cl-fsVps + cl-fsdtVps + + \v^]s - d7s^.PsV<|) . (23) 

2 o 

To achieve Eqs. (|22| ) and (^3]) we have used the lattice sums |^ 

18 

51 = cl^5al3 , (24) 

i=l 

18 ^ 

X] CiaCipCi^CiS = - {5a(35^6 + ^07^/35 + ^aS^p-j) ■ (25) 
i=l 

Eqs. (^) and ( PBj ) are the starting point for the Chapman-Enskog expansion which 
introduces a macroscopic space scale by ri = er and two macroscopic time scales by ti = et 
and t2 = e^t. The ti scale describes fast convection processes whereas on the slower ^2 scale 
the diffusion of vorticity takes place. Thus, the Chapman-Enskog expansion enables us to 
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consider the latter time scales separately. It was for the first time applied to lattice gases 



by Frisch et al. 



We now express the derivatives by means of ri, ti, and 

V = eVi (26) 

dt = edt, + e'dt, , (27) 

and put them into Eqs. (p2|) and (^). If we collect terms of the same order in e we obtain 
on the scale 

dt.Ps = -Vi ■ , (28) 
dtJs = -cl7s{^ips + ZsPs^i^) . (29) 

Eq. (^) is the continuity equation for mass conservation. If one takes the sum over s on 
both sides of Eq. ( PU]) one obtains the "fast" part of the linearized Navier-Stokes equation 
for the total mass current j: 

dtJ=-clY,{Vps + ZsPsVi^) . (30) 

s 

The first term on the right-hand side of this equation is the negative gradient of the pressure 
and the second term the electrostatic body force. 
On the scale we have 

dt.Ps = ^dlp, + dt,V, ■ I + (V?p, + z,Vi • p,V<l) , (31) 

dJs = cl^sdt.ViPs + ^dll + . (32) 



From Eqs. (PE|)-(P2D we see that the transport of the ion densities to equilibrium can be 
either achieved by convection or by diffusion. So, if dtj^js vanishes Eq. (p9|) is solved by the 
Boltzmann distribution, Eq. (^. And by using Eqs. (^) and (^) Eq. (^Tj) for the densities 
simplifies to 

dt.Ps = . (33) 

This equation implies that the fluid is incompressible on the ^2 scale. In this case there is no 
diffusion on the t2 scale because the ionic densities have already come to their equilibrium 
on the faster ti scale. On the other hand, if we assume that the second derivative of ps with 
respect to ti vanishes then Eq. (|3lD becomes a diffusion equation with the diffusion constant 
Dg = c^^'js/^ for ions of type s. The diffusion process relaxes the densities ps {s = +, — ) 
again to the equilibrium distribution (||). 

We still have to discuss Eq. ( |52D for the mass current on the t2 scale. It is reasonable to 
assume that the derivatives of ps and jg with respect to ti are small on the ^2 scale. So we 
may neglect the first two terms in Eq. (P^). Furthermore, we have to sum over s on both 
sides of Eq. (^) in order to obtain the equation of motion for the total current j on the t2 
scale: 

dj = ^ V?j . (34) 
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Moreover, if we combine this equation with Eq. p9|) for variations on the ti scale we obtain 
the following equation: 

dtJ = ^ V^j - cl J2 7s Vp. - cl J2 7s^sPs V<l . (35) 

^ s s 

So we recover Eq. (||) for the total current whereby the kinematic viscosity u is 1/6 and the 
equation for the pressure depends on the parameters 7^, 

P = Csv J2 ^sPs ■ (36) 

s 

Note that for 71 = 72 this is just the equation of state for an ideal gas. 



V. THE SEDIMENTATION VELOCITY 

In this section we present the results for the sedimentation velocity of an array of charged 
spheres in an electrolyte solution. We show that our method gives correct results in that it 
recovers an analytical formula which is valid in the limit of an isolated, weakly charged sphere 
in an unbounded electrolyte. Furthermore, we discuss the dependence of the sedimentation 
velocity on the charge of the spheres. We then compare our results to experimental data. 

Before we show the results for the sedimentation velocity we discuss to what extent the 
calculated electrostatic potentials and the ionic densities around a macroion are influenced 
by lattice artifacts. The latter may have several reasons: We use a simple way to introduce 
the charge of a macroion. We assign its charge Z by distributing the densities on the lattice 
nodes inside the macroion such that the net charge, i.e. the sum over all these nodes, is equal 
to Z. This may have the effect that the effective charge distribution on the surface of the 
macroion is anisotropic because the inner nodes form only an approximative representation 
of a sphere which can, of course, be improved by increasing its radius a. Another drawback 
which is due to the lattice lies in the range of the electrostatic potential. The Debye length, 
measured in lattice units, must be at least larger than one. Otherwise, the simulation 
becomes unstable because one has strong discontinuities in the ionic densities between two 
lattice points close to the surface of the macroion. 

The importance of these artifacts can be checked by determining the radial ionic densities 
Ps(r) around a spherical macroion, i.e. the densities for the ions of type s at a distance r 
from the center of the sphere. In the absence of artifacts Ps{r) should obey 

Ps{r) =Psexp(^-Zs^{r)^ . (37) 

In Fig. ID we show an example of the Ps{r) calculated directly (symbols) and from the right- 
hand side of Eq. (^) (solid lines) for a macroion with a radius a = 4.5 and charge Z = 100. 
The length of the simulation box is L = 80 and the density of the neutral fluid is set to 
Pn = 20. The Debye length is varied by changing the sum of the mean ionic densities 
from J2sPs = 0.0063 {Xd = 5.6) to J2sPs = 0-1 i^D = 1-4) whereby the Bjerrum length 
is /s = 0.4. We can infer from Fig. |I] that the ionic densities as calculated directly agree 
very well with those calculated from the right-hand side of Eq. (pTl) even for a Debye length 
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as small as = 1.4. This means that a = 4.5 is a reasonable choice for the radius of a 
macroion which keeps lattice artifacts small at least for Xd > 1.4. 

Now we consider one charged macrosphere in a unit cell of length L with a fixed posi- 
tion. It represents one particle in an array of spheres on a simple cubic lattice because it 
interacts with its own periodic images. The volume fraction of the macroions is given by 
= (47ra'^)/(3L^). In order to investigate sedimentation phenomena we have to introduce a 
gravitational force in the LB equations. After a time of the order Xg = L'^/ucs [I'es'- effective, 
kinematic viscosity of the fluid in the presence of the spheres) the steady state is reached for 
which one determines the average flow velocity in the unit cell. We divide the latter by the 
average flow velocity of the corresponding neutral system and yield the ratio U/Uq of the 
sedimentation velocities in the charged and the neutral system, respectively. For the larger 
systems considered, we did not wait until the system came to its steady state, because we 
know that the time dependence of the apparent sedimentation velocity t/app(^) is given by 



^app 



W = f/(l-exp(--^)) . (38) 



By computing numerically the time derivative of Eq. (^) one obtains a simple exponential 
function with the two unknown quantities U and Tg which we computed from fits of the 
logarithm of this exponential function. With this procedure it was possible to determine U 
within a time of the order t ^ Ts/20. We checked the accuracy of our fits at = 0.0018 
by comparing them to the exact steady state results, and we obtained identical results for 
U/Uq as a function of na {k = A^^). 

In the following we show that our LB method recovers an analytical result for U/Uq 
which was first derived by Booth [26|, and later slightly modified by Ohshima et al. 



It is valid in the limit of infinite dilution and small charge of the macroions, i.e. a weakly 
charged macroion in an electrolyte with infinite extension. What do we expect in this 
case? Due to the external force the ionic concentrations around the macroions which form 
the electrical double layer deviate from their equilibrium values. The double layer loses its 
spherical symmetry due to the fluid motion which results in an electrical dipole field pointing 
in the direction opposite to the motion of the macroion and thus reduces its sedimentation 
velocity. Booth's calculation starts with the Ansatz 

TJ CO 

- = 1 + E c,Z>^ (39) 

^0 fc=i 

and takes into account only terms in the lowest non- vanishing order in Z'', 

^ = 1 + c^Z' . (40) 

The coefficient C2 can be calculated analytically by solving the electrokinetic equations of 
motion (|l|), (|^), and (|^) whereby the Poisson equation is solved in the Debye-Hiickel limit. 
The final expression for C2 has the following form p6|,|27|]: 



^However, Booth's result agrees with the one of Ohshima et al. for an 1-1-eIectrolyte. We are 
only interested in this special case in this paper. 
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C2 = 



f{Ka) . 



(41) 



72Tia^ri J2s zip, 



Pg denotes the mean density of the ions of type s far away from the center of the macroion. 
7] is the shear viscosity. f{Ka) is a function of exponential integrals of different order n, 



We have determined U/Uq as a function of na for different volume fractions. The radius of 
the macroion is fixed to a = 4.5. Moreover, the diffusivities of the ionic species in the fluid 
are chosen to be Di = 0.165 for the counterions and D2 = 0.25 for the colons. The ratio 
Di/D-2 = 0.66 corresponds to that of D^^/ Dq\ in sodium chloride. The Bjerrum length is 
again set to Ib = 0.4. In order to vary K,a from 0.15 to 1.5 we have to change Z^sPs from 
0.00025 to 0.02211, respectively. This is small compared to the density of the neutral fluid, 
pn = 20. So by changing na we do not change the viscosity of the fluid significantly. 

Fig. ^ shows the results for a surface charge Z = 10 of the macroion. We demonstrate 
below that this value is small enough for the approximation (|28|) to hold. Firstly, we infer 
from Fig. ^ that the relative reduction of the sedimentation velocity due to the charges is 
only of the order of 10~^ at Z = 10 in the yj-range considered. U /Uq exhibits a minimum 
which moves to higher values of na with increasing yj. The occurence of such a minimum is 
reasonable since an increase of na is accompanied with two competing effects: On the one 
hand the electrostatic potential $ becomes stronger due to an increasing salt concentration 
but on the other hand it becomes also more short-ranged because of a decreasing Debye 
length, and thus, it affects only the flow nearby the macroion. The value of the minimum in 
U /Uq decreases with decreasing and seems to move towards the one of the Booth curve for 
(y9 — >• 0. This also holds for the amphtude and the shape of U /Uq for — 0. In order to give 
quantitative evidence that our calculation would recover Booth's result we plot in the inset 
of Fig. H U /Uq as a function of <y9^/^ at na = 0.5, i.e. around the position of the minimum. 
The fit in this figure with a straight line indeed approaches the Booth result, i.e. at = 0. 

Up to now we have shown only the results for U/Uq for a small charge of the macroion 
Z = 10. But it is of course interesting to check up to which values of Z the approximation 
(^) holds. If Eq. (|40|) would be exact, one could renormalize U/Uq as a function of na for a 
given charge Z = Z^m to a new charge Z = Z^^w by multiplying 1 - f//f/o by Z^^^/Z^^^. In 
this way we have renormalized our data for Z = 10 at = 0.00076 to Z = 100 and Z = 130, 
and we compare these data sets in Fig. ^ with the corresponding simulation results for the 
latter two values of Z. We see that we have strong corrections to the results as expected 
from Booth's theory, especially around the minimum in U/Uq. Firstly the amplitude of 
the minimum is underestimated by the renormalized curves and also the position of the 
minimum is at a slightly larger value. To study the corrections to Booth's theory more 
quantitatively we plot in Fig. ^ 1 — U/Uq as a function of Z for = 0.16, 0.5, and 1.5. The 
solid lines in this figure are fits of the form g{Z) = C2Z'^ + c^Z'^ + cqZ^ . From the comparison 
of the different functions g{Z) to the corresponding ones with only the leading term oc Z^ 
(dashed curves in Fig. ^o) we can conclude that the corrections to Booth's theory become 
important for Z > 50. 



E„(a;) = a;"-^ dt t"" exp(-t): 




'e^""^ (3^4 (/«a) - 5£'6(/€a))^ + Se^'^ (^3(^0) - E^i^a)) 



e2«« (4^3(2^0) + 3Ei{2Ka) - 7Es{2Ka)) 



(42) 
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Finally we address the question whether our results are in agreement with a dynamic 
light scattering experiment by Schumacher and van de Ven |2^. They measured the diffusion 
constant D for a system of gold particles in distilled water. Note that D/Dq is equal to U /Uq. 
The radius of the gold particles was around 20 nm and their volume fraction 2 ■ 10~^ %. 
Different data sets were determined by changing the value of na with different salts. We 
consider here the experimental data points measured with sodium chloride which are shown 
in Fig. ^ in comparison to our simulation data at 9? = 0.0046 and at = 0.00076 for Z = 100. 
It is interesting that the experimental data can be very well described by the simulation curve 
for if = 0.0046 although the experiment was done at a very small volume fraction of gold 
particles and the curve for = 0.00076 deviates strongly from the experimental data. More 
systematic experiments, e.g. for different volume fractions, would be necessary to clarify this 
discrepancy. 



VI. CONCLUSIONS 

We have developed a LBM for the simulation of the dynamics of suspensions of charged 
colloidal particles. In this method a set of non-linear, coupled electrokinetic equations is 
solved which consists of convective diffusion equations for the ion densities ps (s = +, — ), 
the linearized Navier-Stokes equations for the mass current j, and the Poisson equation 
for the electrostatic potential Furthermore, a neutral fluid characterized by the density 
Pn can be introduced in order to keep the viscosity in the fluid essentially constant. The 
propagation of ps and j is computed by means of one-particle distribution functions for 
each species s. But in contrast to the normal LBM ps and j are not simply given as zeroth 
and first moments of the nfs, respectively. This is due to the coupling of the ionic part of 
the fluid to the gradient of $ which leads to an additional diffusive term in the propagation 
of the ion densities and a body force term in the propagation of the mass current j. The 
electrostatic potential is determined from the ion densities by means of a Poisson equation 
solver for which we use a successive overrelaxation scheme. Our method is fast, stable and 
easy to implement. We have verified our numerical scheme by means of a Chapman-Enskog 
expansion. 

As an application we applied our method to determine the reduced sedimentation velocity 
U/Uq for an array of charged spheres on a simple cubic lattice. We determined U/Uq as a 
function of the dimensionless parameter Ka for different volume fractions of macroions (p. 
We compared our results with an analytical formula first derived by Booth which is valid in 
the limit of a weakly charged, isolated macroion in an unbounded electrolyte. Booth's theory 
starts with an expansion of U/Uq in powers of the macroion charge Z of which the lowest 
non- vanishing order, Z^ is taken into account. We gave evidence that we recover Booth's 
result in the limit (f 0. For about Z > 50 corrections to the Booth theory become 
important. Up to Z = 130 our simulation data can be well described by an expansion 
up to order Z^. These numerical data could be used to test theories that go beyond the 
Booth level. We mention that very recently also a mode coupling theory with hydrodynamic 
interactions was shown to be in agreement with Booth's theory p9| . This theory is also able 
to consider colloidal systems at finite volume fractions. 

Our LBM for charged colloids is well suited to study their short-time dynamics and 
the flow around macroions. It is not restricted to steady state problems but one can also 
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determine time dependent quantities like the velocity autocorrelation function of a tagged 
macroion in a colloidal suspension. Also particle shapes which are different from a spherical 
one can be easily introduced in our LBM. Moreover, the introduction of walls is rather simple 
which makes it possible to study charged colloidal particles in confined geometries. 
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FIGURES 




5.0 10.0 15.0 20.0 

r 

FIG. 1. The radial ionic densities around a macroion with a = 4.5 for Z = 100 as determined 
directly from Ps{'r) (filled symbols for the counterions and open symbols for the coions) and by 
calculating it from the electrostatic potential (lines) — see text — for the indicated Debye lengths 
Ad. 
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FIG. 2. U/Uo for the volume fractions ip = 0.0046, 0.0018, and 0.00076 as a function of Ka. 

The charge of the macroion is set to Z = 10. The sohd Hne is the result from Booth's theory. The 

inset shows U /Uq as a function of (p^^^ for na = 0.5. The solid line in the inset is the fit function 
0.999269 + 0.00245569?i/3 
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FIG. 3. a) [//C/o at = 0.00076 as a function of kq, for Z = 100 and Z = 130 (filled 
symbols). The open symbols show data for Z = 10 which are renormalized to Z = 100 and 
Z = 130 (see text), b) 1 — U/Uq as a function of Z for the indicated values of Ka. The solid 
lines show the following fit functions: g{Z) = 2.916 • IQ-^^^ _ g 086 • lO'^^Z^ - 1.173 • IQ-^^Z^ 
for Ka = 0.16, g{Z) = 5.182 • lO'^Z^ - 2.810 • lO^^^Z^ + 5.801 • lO'^^Z^ for Ka = 0.5, 
g{Z) = 2.056 • 10^6 - 8.684 • 10"^^ Z*^ + 1.592 • 10"^^ Z^ for Ka = 1.5. The dashed lines show only 
the term oc Z^ of the latter three functions. 
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FIG. 4. U IVq at = 0.00076 and if = 0.0046 as a function of Ka for Z = 100 in comparison to 
experimental data (closed circles). 
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